{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lost Luggage Distribution  Problem\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "In this example, you’ll learn how to use mathematical optimization to solve a vehicle routing problem with time windows, which involves helping a company figure out the minimum number of vans required to deliver pieces of lost or delayed baggage to their rightful owners and determining the optimal assignment of vans to customers.\n",
    "\n",
    "This model is example 27 from the fifth edition of Model Building in Mathematical Programming by H. Paul Williams on pages 287-289 and 343-344.\n",
    "\n",
    "This modeling example is at the advanced level, where we assume that you know Python and the Gurobi Python API and that you have advanced knowledge of building mathematical optimization models. Typically, the objective function and/or constraints of these examples are complex or require advanced features of the Gurobi Python API.\n",
    "\n",
    "**Download the Repository** <br /> \n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Description\n",
    "\n",
    "A small company with six vans has a contract with a number of airlines to pick up lost or delayed baggage, belonging to customers in the London area, from Heathrow airport at 6 p.m. each evening. The contract stipulates that each customer must have their baggage delivered by 8 p.m. The company requires a model to advise them what is the minimum number of vans they need to use and to which customers each van should deliver and in what order. There is no practical capacity limitation on each van. Each van can hold all baggage that needs to be delivered in a two-hour period. To solve this problem, we can formulate an optimization model that minimizes the number of vans that need to be used.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Model Formulation\n",
    "\n",
    "\n",
    "### Sets and Indices\n",
    "\n",
    "$i,j \\in \\text{Locations} \\equiv L=\\{0,1..(n-1)\\}$: Set of locations where $0$ is the index for the single depot -Heathrow airport, and $n$ is the number of locations.\n",
    "\n",
    "$k \\in \\text{Vans} \\equiv  V=\\{0..K-1\\}$: Index and set of vans, where $K$ is the number of vans.\n",
    "\n",
    "$S_k \\in S  $: Tour of van $k$, i.e. subset of locations visited by the van.\n",
    "\n",
    "### Parameters\n",
    "\n",
    "$t_{i,j} \\in \\mathbb{R}^+$: Travel time from location $i$  to location $j$.\n",
    "\n",
    "### Decision Variables\n",
    "\n",
    "$x_{i,j,k} \\in \\{0,1 \\}$: This binary variable is equal 1, if van $k$ visits and goes directly from location $i$ to location $j$, and zero otherwise.\n",
    "\n",
    "$y_{i,k} \\in \\{0,1 \\}$: This binary variable is equal 1, if van $k$ visits location $i$, and zero otherwise.\n",
    "\n",
    "$z_{k} \\in \\{0,1 \\}$: This binary variable is equal 1, if van $k \\in \\{1,2..K\\}$ is used, and zero otherwise.\n",
    "\n",
    "### Objective Function\n",
    "\n",
    "**Number of vans**: Minimize number of vans used.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Minimize} \\quad \\sum_{k = 1}^{K} z_k\n",
    "\\end{equation}\n",
    "\n",
    "### Constraints\n",
    "\n",
    "**Van utilization**: For all locations different from the depot, i.e. $i > 0$, if the location is visited by van $k$, then it is used.\n",
    "\n",
    "\\begin{equation}\n",
    "y_{i,k} \\leq z_{k} \\quad \\forall i \\in L \\setminus \\{0\\}, \\; k \\in V\n",
    "\\end{equation}\n",
    "\n",
    "**Travel time**: No van travels for more than 120 min. Note that we do not consider the travel time to return to the depot.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{i \\in L} \\sum_{j \\in L \\setminus \\{0\\}} t_{i,j} \\cdot x_{i,j,k} \\leq 120 \\quad \\forall k \\in  V\n",
    "\\end{equation}\n",
    "\n",
    "**Visit all customers**:  Each customer location is visited by exactly one van.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{k \\in V}  y_{i,k} = 1 \\quad \\forall i \\in L \\setminus \\{0\\}\n",
    "\\end{equation}\n",
    "\n",
    "**Depot**: Heathrow is visited by every van used. (Note: to improve performance, we diverge from the book by disaggregating this constraint).\n",
    "\n",
    "\\begin{equation}\n",
    "y_{0,k} = z_k \\quad \\forall k \\in V\n",
    "\\end{equation}\n",
    "\n",
    "**Arriving at a location**: If location $j$ is visited by van $k$, then the van is coming from another location $i$.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{i \\in L}  x_{i,j,k} =  y_{j,k}  \\quad \\forall j \\in L, \\; k \\in V\n",
    "\\end{equation}\n",
    "\n",
    "**Leaving a location**: If van $k$ leaves location $j$, then the van is going to another location $i$.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{i \\in L}  x_{j,i,k} = y_{j,k}  \\quad \\forall j \\in L, \\; k \\in V\n",
    "\\end{equation}\n",
    "\n",
    "**Breaking symmetry**: \n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{i \\in L}  y_{i,k} \\geq \\sum_{i \\in L}  y_{i,k+1} \\quad \\forall k \\in  \\{0..K-1\\}\n",
    "\\end{equation}\n",
    "\n",
    "**Subtour elimination**: These constraints ensure that for each van route, there is no cycle. \n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{(i,j) \\in S_k}x_{i,j,k} \\leq |S_k|-1 \\quad \\forall  k \\in K, \\;   S_k \\subseteq L\n",
    "\\end{equation}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Implementation\n",
    "\n",
    "We import the Gurobi Python Module and other Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import math\n",
    "import random\n",
    "from itertools import permutations\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Python 3.7.0 & Gurobi 9.1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input data  \n",
    "We define all the input data for the model. The user defines the number of locations, including the depot, and the number of vans. We randomly determine the coordinates of each location and then calculate the Euclidean distance between each pair of locations. We assume a speed of 60 km/hr, which is 1 km/min. Hence travel time is equal to the distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of locations, including the depot. The index of the depot is 0\n",
    "n = 17\n",
    "locations = [*range(n)]\n",
    "\n",
    "# number of vans\n",
    "K = 6\n",
    "vans = [*range(K)]\n",
    "\n",
    "# Create n random points\n",
    "# Depot is located at (0,0) coordinates\n",
    "random.seed(1)\n",
    "points = [(0, 0)]\n",
    "points += [(random.randint(0, 50), random.randint(0, 50)) for i in range(n-1)]\n",
    "\n",
    "# Dictionary of Euclidean distance between each pair of points\n",
    "# Assume a speed of 60 km/hr, which is 1 km/min. Hence travel time = distance\n",
    "time = {(i, j):\n",
    "        math.sqrt(sum((points[i][k]-points[j][k])**2 for k in range(2)))\n",
    "        for i in locations for j in locations if i != j}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Deployment\n",
    "\n",
    "We create a model and the variables. The decision variables determines the order in which each van visits a subset of custormers, which customer is visited by each van, and if a van is used or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gp.Model('lost_luggage_distribution')\n",
    "\n",
    "# Create variables: \n",
    "\n",
    "# x =1, if van  k  visits and goes directly from location  i  to location  j \n",
    "x = m.addVars(time.keys(), vans, vtype=GRB.BINARY, name='FromToBy')\n",
    "\n",
    "# y = 1, if customer i is visited by van k\n",
    "y = m.addVars(locations, vans, vtype=GRB.BINARY, name='visitBy')\n",
    "\n",
    "# Number of vans used is a decision variable\n",
    "z = m.addVars(vans, vtype=GRB.BINARY, name='used')\n",
    "\n",
    "# Travel time per van\n",
    "t = m.addVars(vans, ub=120, name='travelTime')\n",
    "\n",
    "# Maximum travel time\n",
    "s = m.addVar(name='maxTravelTime')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constraints\n",
    "\n",
    "For all locations different from depot, i.e. $i > 0$, if the location is visited by van $k$, then it is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Van utilization constraint\n",
    "\n",
    "visitCustomer = m.addConstrs((y[i,k] <= z[k]  for k in vans for i in locations if i > 0), name='visitCustomer' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No van travels for more than 120 min. We make a small change from the original H.P. Williams version to introduce a slack variable for the travel time for each van, t[k]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Travel time constraint\n",
    "# Exclude the time to return to the depot\n",
    "\n",
    "travelTime = m.addConstrs((gp.quicksum(time[i,j]*x[i,j,k] for i,j in time.keys() if j > 0) == t[k] for k in vans), \n",
    "                          name='travelTimeConstr' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each customer location is visited by exactly one van"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visit all customers\n",
    "visitAll = m.addConstrs((y.sum(i,'*') == 1 for i in locations if i > 0), name='visitAll' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Heathrow (depot) is visited by every van used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Depot constraint\n",
    "depotConstr = m.addConstrs((y[0, k] == z[k] for k in vans), name='depotConstr' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If location  j  is visited by van  k , then the van is coming from another location  i."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Arriving at a customer location constraint\n",
    "ArriveConstr = m.addConstrs((x.sum('*',j,k) == y[j,k] for j,k in y.keys()), name='ArriveConstr' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " If van  k  leaves location  j , then the van is going to another location  i."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Leaving a customer location constraint\n",
    "LeaveConstr = m.addConstrs((x.sum(j,'*',k) == y[j,k] for j,k in y.keys()), name='LeaveConstr' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Breaking symmetry constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "breakSymm = m.addConstrs((y.sum('*',k-1) >= y.sum('*',k) for k in vans if k>0), name='breakSymm' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Relate the maximum travel time to the travel times of each van"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxTravelTime = m.addConstrs((t[k] <= s for k in vans), name='maxTravelTimeConstr')\n",
    "\n",
    "# Alternately, as a general constraint:\n",
    "# maxTravelTime = m.addConstr(s == gp.max_(t), name='maxTravelTimeConstr')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Objective Function\n",
    "We use two hierarchical objectives:\n",
    "- First, minimize the number of vans used\n",
    "- Then, minimize the maximum of the time limit constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.ModelSense = GRB.MINIMIZE\n",
    "m.setObjectiveN(z.sum(), 0, priority=1, name=\"Number of vans\")\n",
    "m.setObjectiveN(s, 1, priority=0, name=\"Travel time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Callback Definition\n",
    "Subtour constraints prevent a van from visiting a set of destinations without starting or ending at the Heathrow depot. Because there are an exponential number of these constraints, we don't want to add them all to the model. Instead, we use a callback function to find violated subtour constraints and add them to the model as lazy constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Callback - use lazy constraints to eliminate sub-tours\n",
    "def subtourelim(model, where):\n",
    "    if where == GRB.Callback.MIPSOL:\n",
    "        # make a list of edges selected in the solution\n",
    "        vals = model.cbGetSolution(model._x)\n",
    "        selected = gp.tuplelist((i,j) for i, j, k in model._x.keys()\n",
    "                                if vals[i, j, k] > 0.5)\n",
    "        # find the shortest cycle in the selected edge list\n",
    "        tour = subtour(selected)\n",
    "        if len(tour) < n: \n",
    "            for k in vans:\n",
    "                model.cbLazy(gp.quicksum(model._x[i, j, k]\n",
    "                                         for i, j in permutations(tour, 2))\n",
    "                             <= len(tour)-1)\n",
    "\n",
    "\n",
    "# Given a tuplelist of edges, find the shortest subtour not containing depot (0)\n",
    "def subtour(edges):\n",
    "    unvisited = list(range(1, n))\n",
    "    cycle = range(n+1)  # initial length has 1 more city\n",
    "    while unvisited:\n",
    "        thiscycle = []\n",
    "        neighbors = unvisited\n",
    "        while neighbors:\n",
    "            current = neighbors[0]\n",
    "            thiscycle.append(current)\n",
    "            if current != 0:\n",
    "                unvisited.remove(current)\n",
    "            neighbors = [j for i, j in edges.select(current, '*')\n",
    "                         if j == 0 or j in unvisited]\n",
    "        if 0 not in thiscycle and len(cycle) > len(thiscycle):\n",
    "            cycle = thiscycle\n",
    "    return cycle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Verify model formulation\n",
    "\n",
    "m.write('lost_luggage_distribution.lp')\n",
    "\n",
    "# Run optimization engine\n",
    "m._x = x\n",
    "m.Params.LazyConstraints = 1\n",
    "m.optimize(subtourelim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "The optimal route of each van used and the total lost luggage delivery time report follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Print optimal routes\n",
    "for k in vans:\n",
    "    route = gp.tuplelist((i,j) for i,j in time.keys() if x[i,j,k].X > 0.5)\n",
    "    if route:\n",
    "        i = 0\n",
    "        print(f\"Route for van {k}: {i}\", end='')\n",
    "        while True:\n",
    "            i = route.select(i, '*')[0][1]\n",
    "            print(f\" -> {i}\", end='')\n",
    "            if i == 0:\n",
    "                break\n",
    "        print(f\". Travel time: {round(t[k].X,2)} min\")\n",
    "\n",
    "print(f\"Max travel time: {round(s.X,2)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "H. Paul Williams, Model Building in Mathematical Programming, fifth edition.\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
